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Abstract 

Semi-inclusive hadron-production processes are becoming important in high-energy hadron reac- 
tions. They are used for investigating properties of quark-hadron matters in heavy-ion collisions, 
for finding the origin of nucleon spin in polarized lepton-nucleon and nucleon-nucleon reactions, 
and possibly for finding exotic hadrons. In describing the hadron-production cross sections in 
high-energy reactions, fragmentation functions are essential quantities. A fragmentation function 
indicates the probability of producing a hadron from a parton in the leading order of the running 
coupling constant a s . Its Q 2 dependence is described by the standard DGLAP (Dokshitzer- 
Gribov-Lipatov-Altarelli-Parisi) evolution equations, which are often used in theoretical and ex- 
perimental analyses of the fragmentation functions and in calculating semi-inclusive cross sec- 
tions. The DGLAP equations are complicated integro-differential equations, which cannot be 
solved in an analytical method. In this work, a simple method is employed for solving the evo- 
lution equations by using Gauss-Legendre quadrature for evaluating integrals, and a useful code 
is provided for calculating the Q 2 evolution of the fragmentation functions in the leading order 
(LO) and next-to-leading order (NLO) of a s . The renormalization scheme is MS in the NLO 
evolution. Our evolution code is explained for using it in one's studies on the fragmentation 
functions. 
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The DGLAP integrodifferential equations are solved by the Euler's method for the differentiation 
of In Q 2 and the Gauss-Legendre method for the x integral as explained in section |4] 
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This program is used for calculating Q 2 evolution of fragmentation functions in the leading order 
or in the next-to-leading order of a s . Q 2 evolution equations are the timelike DGLAP equations. 
The double precision arithmetic is used. The renormalization scheme is the modified minimal 
subtraction scheme (MS). A user provides initial fragmentation functions as the subroutines 
FFJNI and HQFF in the end of the distributed code FF_DGLAP.f. In FF_DGLAP.f, the subrou- 
tines are give by taking the HKNS07 O) functions as an example of the initial functions. Then, 
the user inputs kinematical parameters in the file setup.ini as explained in section \52\ 
Running time: 

A few seconds on HP DL360G5-DC-X5160. 



1. Introduction 

In the recent years, semi-inclusive hadron-production pro cesses are becoming more and more 
important for studying internal structure of hadrons and heavy-ion reactions. There are three 
ingredients for calculating their cross sections in high-energy reactions with large transverse mo- 
menta for produced hadrons. The first part is on parton distribution functions (PDFs) of initial 
hadrons, the second is on partonic cross sections, and the third is on fragmentation functions 
(FFs) for describing production of hadrons (Qj |2t |3| The perturbative aspect of quantum 

chromodynamics (QCD) has been established for many processes, so that the elementary par- 
tonic cross sections of the second part can be accurately calculated in high-energy processes. 
The PDFs of the first part have been extensively investigated mainly by inclusive deep inelastic 
scattering. Except for extreme kinematical conditions, the unpolarized PDFs are generally well 
determined. For example, one may look at Ref. (0) for the situation of unpolarized PDFs, Refs. 
(0; 0) for polarized PDFs, and Refs. Q fiol) for nuclear PDFs. Because the PDF and partonic 
cross section parts are relatively well known, the only issue is the accuracy of the FFs of the third 
part for calculating precise semi-inclusive cross sections. 

The first estimate for uncertainties of the FFs was done in Ref. S), and its results indicated 
that they have large uncertainties especially for so called disfavored fragmentation functions. 
This fact could add ambiguities to calculated cross sections of high-energy hadron productions 
such as p + p — > n + X and A + A' — > h + X at RHIC and LHC. Here, where p and n indicate 
a polarized proton and a pion, respectively, X indicates a sum over all other hadrons created in 
the reaction, A and A' are nuclei, and h is a produced hadron. Furthermore, the FFs could be 
also used for other studies in searching for exotic hadrons by noting characteristic differences in 
favored and disfavored FFs as pointed out in Ref. 
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These topics suggest that the FFs should be one of most important quantities in describing 
high-energy hadron reactions. There are two important variables in the FFs. One is the energy 
fraction x (or often denoted as z) for a produced hadron from a parton and the other is the hard 
scale Q 2 . Definitions of these quantities are given in Sec. [2] The jc-dependent functions are 
determined mainly from experimental measurements of electron-positron annihilation processes 
e + + e~ — > h + X by a global analysis (Q]; 12; 5 0). On the other hand, Q 2 dependence can 
be calculated in perturbative QCD. The standard equations for describing Q 2 variations are the 
timelike DGLAP evolution equations 

The evolution equations are complicated integro-differential equations, which cannot be solved 
by an analytical method, especially if higher-order corrections are included in the equations. 
They are solved by numerical methods. A popular method is to use the Mellin transformation 
(1121) . in which Q 2 evolution of moments for the FFs is analytically calculated and resulting mo- 
ments are then transformed into x-dependent functions by the inverse Mellin transformation. One 
of other numerical methods is to solve the x-integral part by dividing the x axis into small steps 
for calculating integrals (fl3l [Til) , which is so called brute-force or Euler method. Of course, the 
integral could be calculated by a better method such as the Simpson method or Gauss-Legendre 
quadrature. Another approach is to expand the FFs and DGLAP splitting functions in terms 
of orthogonal polynomials such as the Laguerre polynomials dl5l) . Advantages and disadvan- 
tages of these numerical methods are explained in Ref.(16). There are also recent studies on the 
numerical solution (1171) . 

Although the FFs are important, it is unfortunate that no useful code is available in public for 
calculating Q 2 variations of the FFs. For example, the FFs are calculated in a theoretical model 
at a typical hadron scale of small Q 2 . In order to compare with the FFs obtained by global 
analyses or with experimental data, one needs to calculate the Q 2 evolution. However, one should 
make one's own code or rely on a private communication for obtaining a code since no public 
code is available, although the Q 2 evolution is often used also in theoretical and experimental 
analyses. In this work, we explain our method for solving the DGLAP evolution and a useful 
code is supplied for public use. 

This paper consists of the following. The fragmentation functions and their kinematical 
variables are introduced in Sec. [2] and evolution equations are explained in Sec. [3] Our numerical 
method is described in Sec. |4]for solving the DGLAP Q 2 evolution equations, and a developed 
evolution code is explained in Sec. [5] Numerical results are shown in Sec. [6] by running the 
evolution code, and our studies are summarized in Sec. [7] 



2. Fragmentation functions 

Fragmentation functions are given in the electron-positron annihilation process e + + e~ — > 
h + X, where h indicates a specific hadron. The process is described first by a qq creation by 
<? + e~ — > qq and a subsequent fragmentation, namely a hadron-Zi creation from the primary quark 
or antiquark. The fragmentation function is defined by the hadron-production cross section of 
e + + e- h+X&s ([la): 

fW^J-***"" ->**>, (i) 

o~ fof dx 

where <t, , is the total hadronic cross section. The variable Q 2 is the virtual photon or Z° mo- 
mentum squared in e + e~ — > y (or Z°) and it is expressed by the center-of-mass energy yfs as 
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Q 2 = s. The variable x is the hadron energy Eh scaled to the beam energy yfs/2, and it is defined 
by the fraction: 

E h 2E h 

x = —= — = , (2) 



The fragmentation process is described by the summation of hadron productions from pri- 
mary quarks, antiquarks, and gluons (1191) : 

F\x, Q 2 ) = Yj C * x < ® D * (x > 2 2 )' ( 3 ) 

i 

where C,(x, a s ) is a coefficient function, and it is calculated in perturbative QCD (I20l) . The 
factor a s (Q 2 ) is the running coupling constant, and its expression is given in Appendix A for the 
leading order (LO) and next-to-leading order (NLO). The function Z)' ! (x, Q 2 ) is the fragmentation 
function from a parton i (= u, d, s, ■ ■ •, g) to a hadron h, and it is the probability of producing 
the hadron h, in the LO of a s , from the parton i with the energy fraction x and the momentum 
square scale Q 2 . The notation ® indicates a convolution integral defined by 

f(x)®g(x) = £ jf(y)g{^- (4) 



The fragmentation function is formally given by the expression (12 lh 



r + (o I fc(0, y", X )| /i, x) (/i, x |^(0) 1 



(5) 



where k is the parent quark momentum, the lightcone notation is defined by a* = (a ± a 3 )/ v2, 
the variable x is then given by Jt = plJk + with the hadron momentum ph, and ± is the direction 
perpendicular to the third coordinate. A gauge link needs to be introduced in Eq. (0 so as to 
satisfy the color gauge invariance. It should be, however, noted that a lattice QCD calculation is 
not available for the FFs because the operator-product-expansion method cannot be applied due 
to the fact that a specific hadron h should be observed in the final state with the momentum p/ t . 

An important sum rule of the FFs is on the energy conservation. Since the variable x is 
the energy fraction for the produced hadron, its sum weighted by the fragmentation functions, 
namely the sum of their second moments, should be one. 



dxxD^Q 2 ) = 1. (6) 



The fragmentation function should vanish kinematically at x — 1 and it is expected to be a smooth 
function at small x, so that it is typically parametrized in the form ([ll B 

Dl(x,Ql)=N^(l-xfK (7) 

at fixed Q 2 (= Qfy. Current experimental data are not accurate enough to find much complicated 
x-dependent functional form. In order to calculate the function Z)''(x, Q 2 ) at arbitrary Q 2 , one 
should rely on Q 2 evolution equations and the standard ones are the DGLAP equations in the 
next section. 
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3. Q 2 evolution equations 



The FFs depend on two variables x and Q 2 . The x-dependence is associated with a non- 
perturbative aspect of QCD, so that the only way of calculating it theoretically is to use hadron 
models, because the lattice QCD estimate is not available for the FFs. There are some hadron- 
model calculations by using the expression Eq. (f5]l at a small hadronic Q 2 scale. On the 
other hand, x-dependent functions are determined by global analyses of experimental data mainly 
on e + + e~ — > h +X l[l||2||3i|4|). In the model calculations and also in the global analyses, the Q 2 
dependence or so called scaling violation is calculated in perturbative QCD. 

The Q 2 dependence of the FFs is described by the DGLAP evolution equations in the same 
way with the ones for the PDFs with slight modifications in splitting functions. They are gener- 
ally given by (fTTl [l9l) 



3 r DU X ,Q 2 )= asiQ2) 



ding 2 «i 
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d In Q 



D n Jx,Q') 



2n 



2 , _ <* S {Q 2 ) 



In 



J] P qjVi (x, a.) ® D h q+ {x, Q 2 ) + 2P gq {x, a s ) ® D h g (x, Q 2 ) 
j 

P qg (x, a,) ® D '\^ Q 2) + P ss( x > ">) ® D g(*> G 2 ) 



(8) 



where D q+ (x, Q 2 ) denotes the fragmentation-function combination D h q {x,Q 2 ) + D'~(x,Q 2 ). If 
the sum is taken over the flavor, it becomes the sing let function Dj/je, Q 2 ) = Z q [D'*(x, Q 2 ) + 
D'i(x, Q 2 )], and Nf is the number of quark flavors. The flavor nonsinglet evolution, for example, 
for q - q type function is described by 



a 



d In Q 



■D h q ,{x,Q 2 ) 



In 



Y J Pq 1 q i (x,a s )i 



>D h ^(x,Q 2 ), 



(9) 



where D h q _(x, Q 2 ) 



D h q .{x, Q 2 ) - D h - (x, Q 2 ). The functions P qq (x), P gq {x), P qg (x), and P gg (x) are 
time-like splitting functions, and Pij(x) describes the splitting probability that the parton j splits 
into i with the momentum fraction x. It should be noted that P gq (x) and P qg (x) are interchanged 
in the splitting function matrix for the PDFs. The LO splitting functions are the same as the 
space-like ones; however, there are differences between them in the NLO and higher orders 
(19; 22). Actual expressions of the LO splitting functions are provided in Appendix B. The NLO 



expressions should be found in Ref. ( 19) because they are rather lengthy. 

The DGLAP equations in Eq. (O are coupled integro-differential equations with compli- 
cated x-dependent functions especially if higher-order a s corrections are taken into account. It 
is obvious that they cannot be solved in a simple analytical form. The convolution integral is 
generally expressed by a simple multiplication of Mellin moments, so that the the equations are 
easily solved in the Mellin-moment space. However, the inverse Mellin transformation should 
be calculated by a numerical method in any case to obtain the x-dependent function. Here, we 
solve the DGLAP equations directly in the x space by calculating the x integral in a numerical 
way. Advantages and disadvantages of both methods are discussed in Ref. dig) . 
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4. Numerical method for solving Q 2 evolution equations 



The integro-differential equations of Eqs. © and (|9]l are solved in the following way. Scaling 
violation (Q 2 dependence) of the FFs is roughly given by In Q 2 , which is defined as the variable 
t: 

t = lnQ 2 . (10) 

Because the t dependence is not complicated in the FFs, we do not have to use a sophisticated 
method for solving the differentiation. The following simple method is used for solving the 
differentiation: 

dm _ f{t M )-f(t t ) 

dt At 

Here, the variable t is divided into N, steps with a small interval Af. This method could be called 
Euler method (I23I) . It is also possible to use the Euler method for the integration part by dividing 
the x region into N x steps with the interval Ax (0. However, it is more desirable to use a better 
method since the x dependencies of the FFs and splitting functions are not simple. Here, the 
Gauss-Legendre method is used for calculating the integral over x: 

I g(x) dx =i ° 2^ Wkg(xk), (12) 
Jx o L k=l 

where Xk = [1 + xq + (1 — xo)x' k ]/2 with the zero points xi of the Gauss-Legendre polynomials 
in the region -1 < x' k < +1, Wt are the weights (1241) . and Ncl is the number of Gauss-Legendre 
points. In our previous works (fl3t Il4l) . simpler methods are used for calculating the integral by 
the Euler method and the Simpson's one. Here, we change the method for the Gauss-Legendre 
one for getting more accurate numerical results. 

In the following, only the nonsinglet evolution in Eq. (O is discussed because an extension 
to the general evolution in Eq. © is obvious just by writing down two coupled equations in the 
same way. Substituting Eqs. (fTTT i and (TTZl) into the nonsinglet equation of Eq. (0, we obtain 

D h q - (x m , t M ) = D h q; (x m , t t ) + At Yu 2 Wk J k P W> {xk) D % ' ff ) ■ ( 1 3) 

In previous codes of Q 2 evolution equations in Refs. Sl^t an option is provided to divide 
\nxsj, where Xsj is the Bjorken scaling variable, into equal steps instead of linear-x^y steps 
because the small-x^ region is often important in discussing deep inelastic structure functions. 
However, the small-* part is not as reliable as the PDF case because experimental data do not 
exist at very small x and because of theoretical issues on finite hadron masses and resummation 
effects. The Gauss-Legendre points are taken by considering the linear- x scale at x > 0.1 and 
by the logarithmic-* scale at x < 0.1. If the initial function is supplied at certain Q 2 (= Qq), the 
evolution from t\ - InQ^ to the next point t% = t\ + At is calculated by Eq. (TT~3T >. Repeating this 
step, we finally obtain the evolved FF at fjy+i = In Q 2 . 

The most important and time-consuming part is to calculate the x integrals by the Gauss- 
Legendre quadrature. For the integral from the minimum xq to 1, the splitting functions P q -(xk) 
are first calculated at Ncl points of Xk and they are stored in an array. Then, the fragmentation 
functions are also calculated at x^ and x m , and they are stored in a two-dimensional array. These 
arrays are used for calculating the Gauss-Legendre sum in Eq. (IT~3b . 
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5. How to run the Q 2 evolution code 

We made the numerical evolution code of the FFs by the method discussed in the previous 
section. Its main code (FF_DGLAP.f), a test program (sample. f), and an example of the input file 
(setup.ini) could be obtained upon email request (1251) . There are three major steps for calculating 
the Q 2 evolution of the FFs: 

1. Initial FFs are supplied in the subroutines, FFJNI for gluon (g) and light-quark (u, d, s, u, 
d, s) functions and HQFF for heavy-quark (c, b) functions. 

2. Input parameters for the evolution are supplied in the file setup.ini. These parameters 
are used for calculating two-dimensional (x and Q 2 ) grid data for the FFs in the ranges 
x min < x < 1 and Ql = Q 2 m .„ < Q 2 < Q^ x . 

3. As indicated in the test code (sample. f), the evolved Q 2 value (Q2) and the value of x (X) 
should be supplied for calculating the evolution. The grid data created in the step 2 are 
used for this final step calculation. Therefore, output functions can be obtained at various 
x and Q 2 points without repeating the Q 2 evolution calculations as far as they are within 
the ranges x min < x < 1 and Q 2 } < Q 2 < Q 2 max . 

5.1. Main evolution code 

The main Q 2 evolution code (FFJDGLAP.f) is rather long, so that only the major points are 
explained. First, one needs to supply the initial FFs in the subroutines FFJNI and HQFF, which 
are located in the end of FFJDGLAP.f. The subroutine FFJNI is for gluon (g) and light-quark 
(m, d, s, ii, d, s) functions, and HQFF is for heavy-quark (c, b) functions. As an example, the 
HKNS07 functions ((2) are given. The initial scale for the gluon and light-quark functions is Qq, 
and the scales are the mass-threshold values m 2 and w? for charm and bottom FFs, respectively. 
These scale values are provided in setup.ini, and the initial functions are supplied in analytical 
forms in our main code FFJDGLAP.f. 

Second, input parameters are read from setup.ini, which is explained in Sec. 15.21 They are 
basic parameters: the order of a s , scale parameter of QCD (A), charm and bottom masses m c 
and nib for setting thresholds, number of flavors at the initial scale Q 2 } ; kinematical parameters: 
initial scale Q?, maximum Q 2 value Ql, av and minimum x (x m j„) for making grid data of the 
evolved FFs; parameters to control the numerical integrations: Gauss-Legendre points (Ncl) and 
numbers of In Q 2 = t and x points (N, and N x ). 

Third, the splitting functions are calculated at the points Xk for calculating the summation in 
Eq. ( TOI l. The x^ points are determined by the parameters N x and Ncl- The splitting functions at 
these points are calculated at once in the beginning of this code. In the same way, the initial FFs 
are also calculated at the given points of Xk and x m and they are stored in two-dimensional (k and 
m) arrays. 

Forth, the evolution step of Eq. ( fT3l is repeated in Af, times to obtain the evolved FFs up to 
Q 2 lax from Qq in the range x m - m < x < 1. If Q 2 exceeds the threshold m 2 (or m 2 ), the number of 
flavor is changed accordingly and charm (or bottom) function starts to participate in the evolution 
calculation. During the evolution calculations, two dimensional (x and Q 2 ) grid data are stored 
for calculating the FFs at any point within the ranges x m - m < x < 1 and Qq < Q 2 < Q L maK by 
interpolation. The x and Q 2 values need to be specified in running this main subroutine, and an 
example is proved as a test code (sample. f). 
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5.2. Input file 

The input parameters should be supplied in the file setup.ini for running the main evolution 
routine FFJDGLAP.f, in which the parameter values are read. For example, the following input 
values are used for evolving the HKNS07 FFs in the NLO. Here, the symbol # is for commenting 
out the subsequent line in setup.ini. 



#pQCD ORDER l:LO, 2:NLO 
IORDER= 2 

# DLAM (Scale parameter in QCD) of N f = 4 

# e.g. 0.220 GeV (LO), 0.323 GeV (NLO) in HKNS07 

DLAM= 0.323 # in HKNS07-NLO 

# Heavy-quark mass threshold 

# HQTHRE= m c , m b = 1.43, 4.3 GeV in HKNS07 
HQTHRE= 1.43,4.3 

# Q2 range for making grid files 

# Q02 — > Q2max (note: not the Q2 evolution range) 

Q2= 1.D0, 1558.D+5 # in HKNS07 Library 

# minimum of x 
XMIN= l.D-2 

# NT: # of t-steps for Q 2 dependence 

# NX: # of x-steps for xD(x) 

# NGLI: # of steps for Gauss-Legendre integral 
NT= 580 

NX= 160 
NGLI= 32 

# NF at the initial scale Q02 

NF= 3 (14) 

The parameter IORDER indicates the LO or NLO of a s . Both LO and NLO evolutions are 
possible so that the order of a s should be IORDER=l or 2. The scale parameter A should be 
supplied in the case of four flavors (A4). It is then converted to the three and five flavor values 
(A3 and A5) within the evolution code (26). The charm and bottom functions appear as finite 
distributions above the threshold values Q 2 > m 2 (or nth. The HQTHRE values are these heavy- 
quark thresholds (1271) . 

The kinematical regions of x and Q 2 are specified by the parameters, Qq (=Q 2 ni „), Qm ax > 
and x m ,„. The Q 2 is the initial Q 2 scale, and Q 2 1WX is the maximum Q 2 for calculating the FFs. 
Any Q 2 values can be chosen as long as perturbative QCD calculations are valid, which means 
that small £>o an d Qmax values are not favorable particularly in the region Q 2 < 1 GeV 2 where 
pQCD calculations do not converge easily due to the large running coupling constant a s . The Q 2 
maximum Q 2 max = 1.558 x 10 8 GeV 2 is used in making the HKNS07 library for their FFs (0). It 
is chosen so that two grid points are close to the charm and bottom threshold values m 2 and m?. 
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If one needs to calculate the evolution only to 2 2 -100 GeV 2 , one does not have to take such a 
large Q 2 mx , and 2„ aA = 100 GeV 2 is enough. We also should note that a very small value of jc m „, 
is not favored because some FFs become negative, which is not physically allowed in principle. 
This occurs due to a singular behavior of a time-like splitting function. In order to cure this issue, 



more detailed studies are needed by including resummation effects (1281) . 

The parameter NT is the number of points of In Q 2 in the range In Q 2 min < In Q 2 < In Ql, ax , 
and NX is the number of points of x in the range x m - m < x < 1 . The NGLI is the Gauss-Legendre 
points for numerical integration. In the example of Eq. ([14}, the number of t = In Q 2 points is 
580, the one of x is 160, and the one of the Gauss-Legendre points is 32. They are selected by 
looking at evolution results by varying their values. Such studies are explained in details in Sec. 
[6] The NF is the number of flavors at Q 2 r and it is usually three as given in Eq. (11411 . 

5.3. Sample code 

A test code (sample. f) is provided as an example for running the main evolution code. The 
main evolution routine GETFF(Q2,X,FF) is called by supplying Q 2 and x values. Evolved FFs 
are returned to FF(I), (I=-5, 5): 

FF(-5) = D\{x, Q\ FF(-4) = D^x, Q 2 ), FF(-3) = D^x, Q 2 ), FF{-2) = D h u {x, Q 2 ), 
FF(-\) = D'p,Q 2 ), FF(Q) = D h g (x,Q\ FF(1) = D h d (x, Q 2 ), FF(2) = D h u {x, Q 2 ), 



FFO) = D h s {x,Q 2 ), FF{A) = D h c (x,Q 2 ), FF(5) = D h b (x, Q 2 ). 



(15) 



6. Results 

Q 2 evolution results of FFs are shown in Fig. [TJby taking the initial functions of n + as the 
HKNS07 (Hirai, Kumano, Nagai, Sudoh) parametrization in 2007 (2). The initial functions are 
provided at Q 2 -\ GeV 2 for g, u, d, and s FFs and for c and b at Q 2 — m 2 and m?, respectively. 
The evolution has been calculated in the NLO and with the scale parameter Aqcd =0.323 GeV 
in the running coupling constant. The used numbers of steps are N,=560, N X -16Q, and Ngl=32 
for calculating the evolutions. In Fig. [2] the Q 2 evolution results are shown as a function of Q 2 
at fixed x points (x = 0.1 and 0.4). The same input parameters are used in setup.ini, which was 
used in obtaining the results in Fig. Q~| for running the code. 

The input file setup.ini for calculating the evolution of the NLO HKNS07 functions is given 
in Eq. 0141 . The light-parton (g, u, d, s, u, d, s) FFs are supplied at the initial scale Q 2 y which is 
assigned to be the Q minimum Q 2 min . The evolved Q 2 value (10, 100, or 10000 GeV 2 ) needs to 
be supplied when running the code, for example, sample. f. 

Next, evolution results are shown by varying the parameters Nql, N x , and A^,, which affect 
the evolution accuracy. First, the Q 2 evolution results are calculated at Q 2 -\QQ GeV 2 by using 
the HKNS07 parametrization for the initial functions and the parameters ^^ = 100, A^ v =500, 
and Af,=500. Then, they are considered to be "standard" functions in showing ratios with other 
evolution results. In the input setup.ini file, Q 2 mx -\QQ GeV 2 is taken because the larger-Q 2 
region is not necessary. 

First, Ngl is varied as 10, 20, and 40 in order to find its dependence on evolution results in 
checking evolution accuracy. The evolved functions are then used for calculating ratios with the 
standard evolution by Df(x,Q 2 = 100GeV 2 ) NGL , N - 5 oo,N,=50o/Df(x,Q 2 = 100GeV 2 ) Woi= 
The ratios are shown in Fig. [3] for the fragmentation functions of g, u, d, c, and b. The quark 
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Figure 1: Q 2 evolution of HKNS07 fragmentation functions. The initial g, u, and d FFs are supplied at the scale Qq=1 
GeV 2 , and the c and b functions are at m 2 and b 2 b . They are evolved to the scale Q 2 =10, 100, and 10000 GeV 2 by the 
time-like DGLAP evolution equations in the NLO (MS ) by using the code developed in this work. The explicit parameter 
values are listed in Eq. )14t . 
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Figure 2: Q 2 dependence of the fragmentation functions at flexed x (= 0.1 and 0.4) . The initial g, u, and d FFs are 
supplied at the scale 2q=1 GeV 2 , and the c and b functions are at m 2 and b 2 . They are evolved to the scale Q 2 =10000 
GeV 2 by the time-like DGLAP evolution equations in the NLO (MS ) by using the code developed in this work. The 
explicit parameter values are listed in Eq. dl4l . 



functions are evolved accurately except for the region close to x — 1 even with a small number 
of Gauss-Legendre points such as Ncl = 10. However, the gluon evolution depends much on the 
choice of Ncl, and the results indicate that Ncl ^ 20 needs to be taken for getting the evolution 
accuracy better than about 0.3%. This is the reason why NGLI = 32 is used in calculating the 
evolutions in Fig. [TJ 

Second, the dependence on N x is shown by fixing the other parameters at A^gl=100 and 
N t -5Q)Q. The evolved functions are used for taking ratios with the standard evolution results with 
N x =500, and N r -500 in the same way with Fig. [3] The input parameter N x is the 
number of points in x from x„„„ to one. For example, if jc m ,„=0.01 and N x =500 are taken, 250 
points are given for the logarithmic x in the region 0.01 < x < 0.1 and another 250 points for 
the linear x in 0.1 < x < 1. If ;t m ,„=0.001 and N X -6QQ are taken, we have 400 (200) points in 
0.001 < x < 0.1 (0.1 < x < 1). It is changed as N x -20, 50, and 200 to show variations in the 
evolved functions, and results are shown in Fig. |4] In general, there are large differences in the 
large-x region in all the FFs. In particular, the evolved functions are not reliable at x > 0.7 if 
N X =2Q is taken. As the number increases as N X -5Q and 200, they become reliable except for the 
extremely large-jc region (x > 0.9). From these studies, Af t =160 is taken, for example, in Fig. Q] 
as a reasonable choice. 

Third, we show A^ f dependence in Fig. [5] It is varied as Af ; =100, 200, and 300 by fixing other 
parameters at Ncl = 100 and N x - 500. If A^, is small, evolved distributions are not accurate at 
large x, especially in the gluon fragmentation function. A large number of points should be taken 
for Af, for getting a converging function within a few percent level of accuracy, and A^,=580 is 
taken in Fig. [1] However, if Qj naK is small, accurate evolution results can be obtained by taking 
smaller N t . 

A typical running time for obtaining the evolutions in Fig. [TJis 4 seconds by using g95 on 
the CPU (Dual-Core Intel Xeon 2.66 GHz) with Mac-OSX-10.5.8, so that the code is efficient 
enough to be used on any machines for one's studies on the fragmentation functions. 
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Figure 3: Evolved fragmentation-function ratios (x)n GLi n x =5Q0. w,=50o/-D J 1 (*)jv oi =ioo, JV,-=500, iV,=500 are shown for 
N G l=W, 20, and 50 at g 2 =100 GeV 2 . The initial functions are the HKNS07 ones at g 2 =l GeV 2 for g, u, and d, at 



Q 1 = m 2 for c, and at Q 2 = mjj for b. 
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Figure 4: Evolved fragmentation-function ratios (-c)(v Gt =loo, jv,. jV,=50o/f f WA' ci =l00, jV v =500, jv,=soo are shown for 
N x =20, 50, and 200 at g 2 =100 GeV 2 . The other conditions are the same as the ones in Fig. [3] 
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gure 5: Evolved fragmentation-function ratios D? (*)iv G1= ioo, n x =500, Ntl^t ( x )n gl =ioo, #,=500, Af,=500 are shown for 
v =100, 200, and 300 at Q 2 =100 GeV 2 . The other conditions are the same as the ones in Fig. [3] 
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7. Summary 



The fragmentation functions are used in describing hadron-production cross sections at high 
energies. The FFs are described by two variables x and Q 2 . The Q 2 dependence of the FFs is 
calculated in perturbative QCD and they are described by the DGLAP evolution equations. In this 
work, the Q 2 evolution equations are numerically solved and a useful evolution code is provided 
so that other researchers could use it for their own studies. The variables x and In Q 2 are divided 
into small steps, and the evolution is numerically calculated by using the Euler method and the 
Gauss-Legendre quadrature. We showed that the evolution is accurately calculated except for 
the extremely large-x region by taking reasonably large numbers of the Gauss-Legendre points 
(Nql), x steps (N x ), and t — In Q 2 steps (N t ). Our evolution code can be obtained upon request 



d25l) for using one's studies on the Q 2 evolution of the FFs 
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Appendix A. Running coupling constant 

The running coupling constants in the leading order (LO) and next-to-leading order (NLO) 



are 



Att 

a ^ Ql) -^J^y (16) 



A,ln(e 2 /A 2 ) 



A lnln(e 2 /A 2 ) 



1 - 



$ln(G 2 /A 2 ) 



(17) 



where A is the QCD scale parameter, and /3o and f3\ are given by 

11 4 34 , 10 

A) = —C G - -T R N f , = —C 2 C - —C c Nf - 2C F N f , (18) 

with the color constants 

N 2 - 1 1 
C A = N C , C F = ^-, T R = ~. (19) 

In the NLO, MS is used for the renormalization scheme. 

Appendix B. Splitting functions 

The splitting functions are expanded in a s : 

Pij(x, or,) = Pf/ix) + ^lp^(x) + ■■■, (20) 
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where P®>(x) and P<ff(x) are LO and NLO splitting functions, respectively. Splitting functions 

J L 

in the LO are the same as the ones for describing the PDF evolution ( 13): 

1 + x 2 3 



(l-x) + 2 
T R [ x 2 + (1 - xf ] , 
1 +(1 - xf 



+ ~ 6(1 -x) 



C F - 

2C C 



x 
x 



1 



(l-*) 4 



+ x(\ - x) + 



11 

L2 



1 A^Te 
3 C G 



6(1 -x) 



(21) 
(22) 
(23) 
(24) 



The only point one should note is that the splitting functions P qg and P gq are interchanged in the 
matrix of Eq. ([8} from the PDF evolution. However, the spacelike and timelike splitting functions 
for the PDFs and FFs, respectively, are different in higher-order of a s as shown in Refs. 
The quark-quark splitting function in the NLO is given by 



pQ) _ p (i) , p (i) 

q*q* ~ mi q&) 



x..(pV0) 4. p v m\ + psm , pS(D 

°']\ r qq + r qq > + r aa + ' n 



qq 



qq 



(25) 



where the functions P qq V) and are given in Ref. the function P^ 1 ' (= P 3 ^) can be 



derived from the relation P'„ - , 
and they are provided in Sec. 6.1 of Ref. fll9| 



qq 

P V JP + P V ,^ + N/ipSJp + P S qq 1} ). These expressions are lengthy 
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sample.f 

c 

PROGRAM SAMPLE ! 2011-11-24 
c 

C X DEPENDENCE OF FFs 

IMPLICIT REAL*8(A-H,0-Z) 
PARAMETER (NSTEP=200) 
CHARACTER* 1 Q2PR0G_END 
CHARACTER*9 Q2_0UT 
CHARACTER*22 FNAME 
DIMENSION FF(-5:5) ,Z(600) 

COMMON /XMIN/XMIN ! Setting in the setup.ini 

CALL FF_DGLAP() ! Making arrary of FF(XMIN: 1 .DO, Q2_ini : Q2_max) 

200 WRITE(*,fmt=' (a) ') "Q~2= "; READ ( * , * ) Q2 
WRITE (Q2_0UT, ' (1PE9.3) ') Q2 
FNAME= ' Q2= ' //Q2_0UT// ' _GeV2 . dat ' 
OPEN (unit=23 , f ile=FNAME , FORM= ' formatted ' ) 

DLMIN=DL0G10 (XMIN) 

ZLSTEP= (DL0G10 ( 1 . DO) -DLMIN) /DFLOAT (NSTEP) 
DO I=1,NSTEP+1 

DLOGZ=DFLOAT ( I- 1 ) *ZLSTEP+DLMIN 

Z(I)=10.D0**(DL0GZ) 
END DO 

C FOR pi~+, FF(I), 1= 0:g, l:d, 2:u, 3:s, 4:c, 5:b 
DO 1=1, NSTEP 



CALL GETFF(Q2,Z(I) ,FF) ! Getting FF(Z,q~2) 





WRITE(23,1010) Z(I), Z(I)*FF(0), 


gluon 


+ 


Z(I)*FF(2) , 


up 


+ 


Z(I)*FF(1), 


down 


+ 


Z(I)*FF(3), 


strange 


+ 


Z(I)*FF(4) , 


charm 


+ 


Z(I)*FF(5) 


bottom 



END DO 

WRITE(*,fmt=' (a) ') 
+ "Do you finish the FF Q2 evolution ? (y/n) " 

READ(*,*) Q2PR0G_END 

IF((Q2PR0G_END(1:1) .EQ . >n') .OR. (Q2PR0G_END(1 : 1) .EQ . 'N')) GOTO 200 

1010 F0RMAT(1X,9(1PE16.7)) 
CL0SE(23) 
END 

c 
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TEST RUN OUTPUT 

Running the distributed sample code (sample. f) together with the main Q 2 evolution sub- 
routine (FFJDGLAP.f) and the input file (setup.ini), we obtain the following output for Q 2 -100 
GeV 2 . The following functions corresponds to the curves at Q 2 -100 GeV 2 in Fig. Q] 



X 




xJr t 


xTT u 




xD". 

a 




xir s 




xD* 




X b 


1.0000000E- 


■02 


1.7387272E+00 


1.5206099E+00 


7 


1337454E- 


■01 


7.1374282E-01 


1 


.2897181E+00 


2 


.5241504E+00 


1.0232930E- 


■02 


1.7469241E+00 


1.5297365E+00 


7 


2321764E- 


■01 


7.2358058E-01 


1 


.2947968E+00 


2 


.5089220E+00 


1.0471285E- 


■02 


1.7544533E+00 


1.5384281E+0O 


7 


3265629E- 


■01 


7.3301393E-01 


1 


.2995278E+00 


2 


.4936096E+00 


1.0715193E- 


■02 


1.7613271E+00 


1.5466907E+00 


7 


4169702E- 


■01 


7.4204939E-01 


1 


.3039165E+00 


2 


.4782151E+00 


1.0964782E- 


■02 


1.7675595E+00 


1. 55453 11E+00 


7 


5034733E- 


■01 


7.5069444E-01 


1 


.3079690E+00 


2 


.4627404E+00 


1.1220185E- 


■02 


1.7731 63 1E+00 


1.5619553E+00 


7 


5861384E- 


■01 


7.5895574E-01 


1 


.3116905E+00 


2 


.4471874E+00 


1.1481536E- 


■02 


1.7781510E+00 


1.5689699E+00 


7 


6650367E- 


■01 


7.6684037E-01 


1 


.3150869E+00 


2 


.4315579E+00 


1.1748976E- 


■02 


1.7825365E+00 


1.5755812E+00 


7 


7402396E- 


■01 


7.7435550E-01 


1 


.3181639E+00 


2 


.4158538E+00 


1.2022644E- 


■02 


1.7863321E+00 


1.5817955E+00 


7 


8118145E- 


■01 


7.8150786E-01 


1 


.3209270E+00 


2 


4000771E+00 


1.2302688E- 


■02 


1.7895498E+00 


1.5876182E+00 


7 


8798236E- 


■01 


7.8830368E-01 


1 


.3233809E+00 


2 


.3842291E+00 


8.9125094E- 


■01 


2.7586761E-05 


7.8509439E-03 


1 


0318809E- 


■06 


1.0318809E-06 


4 


.27834 16E-05 


1 


.7161909E-06 


9.1201084E- 


■01 


1.4513488E-05 


5.0877079E-03 


2 


3777347E- 


■07 


2.3777347E-07 


1 


.6972984E-05 


4 


.98288 19E-07 


9.3325430E- 


■01 


6.3142094E-06 


2.8910396E-03 


3. 


2396330E- 


■08 


3.2396330E-08 


5 


.0850503E-06 


9 


.9610221E-08 


9.5499259E- 


■01 


1.9415364E-06 


1.2925565E-03 


7 


2014456E- 


■10 


7.2014458E-10 


9 


.1254976E-07 


1 


.0039050E-08 


9.7723722E- 


■01 


2.5705072E-07 


3.2218983E-04 




2.2238409E-10 


-2.2238409E-10 


4 


.7361837E-08 


1 


.7555729E-10 
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